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Abstract 

Explicit  volume  averaging  procedures  are  used  to  motivate  a  gradient-type  description  of  single  crys¬ 
talline  elastoviscoplasticity.  Upon  regarding  local  elastic  and  plastic  deformation  gradients  within  the 
crystal  as  continuously  differentiable  fields,  we  arrive  at  a  three-term  multiplicative  decomposition  for  the 
volume-averaged  deformation  gradient,  consisting  of  a  recoverable  elastic  term  associated  with  the  average 
applied  stress  and  average  lattice  rotation,  an  inelastic  term  associated  with  the  average  plastic  velocity 
gradient,  and  a  (new)  third  term  reflecting  the  presence  of  the  residual  microelastic  deformation  gradient 
within  the  volume  and  providing  a  representation  of  the  kinematics  of  grain  subdivision  via  formation  of 
low-angle  subgrain  boundaries,  for  example.  A  variant  of  the  classical  Eshelby  stress  tensor  provides  the 
driving  force  for  homogenized  viscoplastic  flow,  with  slip  resistances  dictated  by  densities  of  geometrically 
necessary  and  statistically  stored  dislocations.  Distinctive  features  of  the  continuum  model  include  coupling 
of  internal  elastic  strain  energy  densities  associated  with  residual  and  applied  stresses,  dependency  of  the 
single  crystalline  effective  elastic  moduli  upon  evolution  of  lattice  substructure,  and  a  characteristic  length 
potentially  based  upon  both  the  size  of  the  crystal  element  used  in  volume  averaging  and  the  grain  sub¬ 
division  measure. 
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1.  Introduction 

A  rather  large  number  of  single  crystal  and  polycrystal  plasticity  models  have  been  proposed 
incorporating  higher  than  first-order  gradients  (e.g.,  strain  gradients,  lattice  curvature,  gradient- 
based  dislocation  density  measures)  of  deformation  in  the  material  response  [1-6].  The  higher 
order  gradients  render  these  classes  of  models  non-local,  effectively  injecting  a  length  scale  (e.g., 
normalization  constant)  into  the  formulation  on  dimensional  grounds,  the  value  of  which  is 
presumably  associated  with  a  characteristic  dimension  of  the  microstructure,  such  as  the  evolving 
dislocation  cell  size  in  ductile  FCC  crystals,  for  example  (cf.  [7]). 

Higher  order  deformation  gradients  have  been  included  within  plasticity  theories  for  a  variety 
of  different  reasons;  we  elaborate  briefly  on  a  few  of  these  here.  Their  inclusion  has  permitted 
resolution  of  numerical  difficulties  associated  with  the  solution  of  boundary  value  problems  of 
strain  softening  materials  in  which  strain  localization  occurs  (cf.  [8,9]).  Gradient  based  approaches 
have  also  been  used  to  model  dislocation  dynamics  and  pattern  formation  [10-12]  and  to  describe 
single  and  periodic  shear  bands  in  single  crystals  (cf.  [13])  and  poly  crystals  [12].  Other  recent 
applications  of  non-local  theories  include  characterization  of  stress  and  strain  fields,  without 
singularities  in  the  field  variables,  at  dislocation  cores  and  crack  tips  [14,15]  and  modeling  the 
evolution  of  the  plastic  spin  tensor  in  macroscopic  finite  deformation  plasticity  theory  [16,17]. 
Perhaps  the  most  recurrently  reported  motivation  for  use  of  gradient  theories  of  plasticity  has 
been  representation  of  the  observed  trend  of  increasing  strength  with  decreasing  size  of  considered 
volume  or  micro  structural  features.  Often  cited  is  the  Hall-Petch  relation,  in  which  hardness 
properties  (i.e.,  yield  stress  and  cleavage  strength)  increase  with  decreasing  grain  size  in  poly¬ 
crystals  (specifically  an  inverse  square-root  dependence  [18,19]),  a  phenomenon  that  classical  local 
plasticity  theory,  being  devoid  of  a  material  length  scale,  is  unable  to  capture.  Shu  and  Fleck  [20] 
and  Forest  et  al.  [21]  used  couple  stress  theories  (pioneered  by  Cosserat  and  Cosserat  [22,23])  to 
characterize  Hall-Petch  behavior  in  bicrystals  and  polycrystals,  respectively.  Fleck  et  al.  [24] 
employed  a  couple  stress  model  of  strain  gradient  plasticity  to  describe  an  increase  in  flow  stress 
with  decreasing  diameter  of  twisted  thin  copper  wires.  Shu  and  Fleck  [25]  and  Hwang  et  al.  [26] 
used  variations  of  the  same  strain  gradient-couple  stress  theory  [27]  to  capture  an  observed  in¬ 
crease  in  hardness  with  decreasing  indentor  size  in  pure  metals. 

Many  theories  have  included  a  higher  order  gradient  of  elastic,  plastic,  or  total  deformation 
in  the  expression  for  the  yield  stress,  slip  system  hardening  rate,  flow  stress,  or  the  backstress, 
with  details  of  incorporation  of  gradients  in  the  material  response  functions  varying  widely 
among  different  models.  In  many  cases  the  higher  order  gradients  are  associated,  upon  invo¬ 
cation  of  differential-geometric  arguments,  with  the  density  of  “geometrically  necessary  dislo¬ 
cations”  (GNDs)  in  single  crystals  (cf.  [28-37])  and  in  polycrystals  (cf.  [38,16,17,9,39]).  These 
dislocations  are  required  to  sustain  the  compatibility  of  the  total  deformation  gradient,  in  con¬ 
trast  to  the  “statistically  stored”  dislocations  (SSDs)  [40]  that  accumulate  under  homogeneous 
plastic  flow. 

Physical  experiments  have  demonstrated  how  the  processes  of  grain  subdivision  and  dis¬ 
location  substructure  formation  substantially  influence  slip  system  activity,  strain  hardening, 
stored  lattice  energy,  and  texture  evolution  in  single  and  polycrystals  [41-47].  Local  models 
have  been  proposed  that  explicitly  embed  subdivision  and  related  dislocation  substructure 
effects  into  the  single  crystal  kinematics  [47-50]  and  the  hardening  and  grain  interaction  laws 
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of  polycrystal  plasticity  theory  [51,52],  without  explicitly  considering  the  higher  order  gradi¬ 
ents  of  deformation  associated  with  geometrically  necessary  defects.  Non-local  solutions  in  terms 
of  Green’s  functions  [14,53]  have  also  been  developed  to  address  defect  substructure  within  grains 
and  its  influence  on  the  micro-stress  fields  and  commensurate  hardening  behavior.  Also 
noteworthy  is  the  recent  approach  of  Ortiz  and  co-workers  [7,54]  for  modeling  additional 
internal  degrees  of  freedom  associated  with  dislocation  substructure  via  a  non-local  sequential 
lamination  theory.  This  approach  permits  dislocation  substructure  development,  when  ener¬ 
getically  favorable,  to  occur  within  single  crystals  even  under  uniform  monotonic  macroscopic 
deformation,  in  contrast  to  many  gradient-based  approaches  requiring  heterogeneity  of  the 
deformation  gradient  field  to  drive  evolution  of  GNDs.  Such  lamellar  dislocation  structures 
have  been  observed  experimentally  within  ductile  single  crystals  under  monotonic  loading  at 
large  deformations  [42,43,55]  and  are  thought  to  influence  the  flow  stress  across  a  wide  range 
of  temperatures  and  strain  rates  [44,56].  Carstensen  et  al.  [57]  viewed  evolution  of  plastic  flow 
from  the  standpoint  of  a  constrained  energy  minimization  problem  and  remark  how  heteroge¬ 
neous  plastic  flow  (e.g.,  subdivision  or  localization)  may  result  from  a  lack  of  convexity  of  the  free 
energy. 

Also  observed  within  pure  ductile  metals  and  certain  alloys  at  large  deformations  and/or  high 
temperatures  are  long  range  internal  stress  fields  associated  with  misoriented  subgrain  bound¬ 
aries  [58-60],  These  internal  stress  fields — attributed  to  misorientations  that  develop  among 
neighboring  subgrains — occur  with  a  periodicity  on  the  order  of  the  subgrain  size  (cf.  [43]).  The 
study  of  Gibeling  and  Nix  [59]  found  that  while  local  internal  stress  fields  in  the  unloaded 
configuration  (i.e.,  residual  stresses)  are  typically  quite  small  compared  to  the  average  applied 
stress,  the  applied  loads  can  alter  the  arrangement  of  subgrain  walls  such  that  the  stress  fields  of 
adjacent  dislocations  do  not  cancel,  thereby  biasing  the  internal  stress  fields  to  a  substantial 
degree. 

The  remainder  of  this  work  is  organized  as  follows.  Kinematics  and  balance  laws  are  discussed 
in  Sections  2  and  3,  respectively,  from  the  perspective  of  multiscale  volume  averaging.  Section  4 
more  fully  develops  constitutive  relations  within  the  context  of  a  higher  order  gradient,  single 
crystal  plasticity  theory.  In  the  interest  of  brevity,  thermal  effects  (i.e.,  temperature  rates  and  heat 
fluxes)  and  dynamic  effects  (i.e.,  acceleration  and  body  forces)  are  often  neglected. 

We  employ  the  following  notation.  Vector  and  tensor  quantities  are  typically  represented  with 
boldface  type,  while  scalars  and  individual  components  of  vectors  and  tensors  are  written  in 
italics.  The  index  notation  is  often  used  for  clarity,  following  the  Einstein  summation  convention 
and  distinguishing  between  covariant  (subscript)  and  contravariant  (superscript)  components. 
Current  configuration  indices  are  written  in  lower  case  Latin,  reference  configuration  indices  in 
upper  case  Latin,  and  intermediate  configuration  indices  are  written  using  Greek  symbols.  Jux¬ 
taposition  implies  summation  over  two  repeated  adjacent  indices  (e.g.,  (AB)fl6  =  AacBch).  The  dot 
(scalar)  product  of  vectors  is  represented  by  the  symbol  (e.g.,  a  •  b  =  aagahbh,  with  gab  com¬ 
ponents  of  the  metric  tensor).  Angled  brackets  denote  a  dual  (scalar)  product  (e.g.,  for  second- 
rank  tensors,  (A,  B)  =  tr(AB)  =  AahBba,  and  for  contra  covector  pairs,  (a,  b)  =  a aba).  The  colon 
denotes  contraction  over  repeated  pairs  of  indices  (e.g.,  A  :  B  =  tr(A7 B)  =  AahBah  and 
C  :  E  =  CahcdEcd).  The  symbol  “<g>”  represents  the  tensor  (outer)  product  (e.g.,  (a  <g>  b)"/’  =  aabh). 
Superposed  -1,  T,  and  denote  inverse,  transpose,  and  material  time  derivative  operations, 
respectively.  Additional  notation  is  clarified  as  needed  later  in  the  text. 
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2.  Multiscale  kinematics 

A  fundamental  difference  between  the  model  developed  in  the  present  work  and  the  numerous 
strain  gradient  approaches  already  cited  is  our  usage  of  rigorous  volume  averaging  procedures  to 
characterize  the  kinematics  of  elasto(visco)plastic  deformation.  Continuum  elements  (i.e.,  do¬ 
mains  for  volume  integration)  representing  a  deforming  crystal  or  region  within  are  shown  in  Fig. 
1.  The  symbols  vie{,  vmt,  and  rcur  denote  volume  elements  in  reference,  intermediate,  and  current 
configurations.  Dimensions  of  /'ref  are  assumed  to  adhere  to 

"C  fref  ^  La,  (1 ) 

with  uq  and  La  the  lattice  parameter  and  average  grain  diameter,  respectively,  and  with 
fref  =  We  also  frequently  invoke  the  following  notation  for  configurations  at  multiple  length 
scales:  crer  =  bre f  c  B[e f,  f]nt  =  bmt  C  BmX,  and  rcur  =  bcm  c  Bcu r,  where  global  reference,  interme¬ 
diate,  and  current  configurations  of  the  entire  macroscopic  body  are  labeled  BIc{,  BmU  and  Bcm.  All 
configurations  will  be  described  in  more  detail  later  in  the  text.  Local,  scalar  “differential”  volume 
elements  within  reference,  intermediate,  and  current  configurations  are  labeled  drref,  dtv,  and 
df'cur-  respectively.  Differential  elements  are  required  to  be  smaller  than  their  associated  volume 
elements,  but  not  necessarily  infinitesimal  in  size. 

The  motion  of  points  of  the  body  from  the  reference  to  current  configuration  is  written 

xa  =  <pa{XA,t)  (a,  A  =  1, 2, 3),  (2) 

where  xa,  XA,  and  /  denote,  respectively,  current  configuration  coordinates,  reference  configura¬ 
tion  coordinates,  and  time.  The  local  motion  tpa  is  assumed  to  be  continuously  differentiable 
within  rcur.  Local  intermediate  configuration  coordinates  are  also  assumed  to  be  continuous  within 
v.  These  are  denoted  by  xf  =xx(XA,t)( a  =  1,2,3).  Further  remarks  on  the  availability  of  con¬ 
tinuous,  single-valued  coordinates  xf  follow  later. 

Basis  vectors  in  each  configuration  are  labeled  as 


Fig.  1.  Configurations,  differential  volumes,  and  coordinate  systems  for  crystal  volume  element:  (a)  reference,  (b)  inter¬ 
mediate,  and  (c)  current. 
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G"  '  cXA  ’  87  ~  Qx  J  ’  8,1  ~  dxa  ' 

Basis  covectors  (i.e.,  dual  vectors)  are  introduced  as 

GA  =  dXA,  g7  =  dF,  g"  =  dxa, 
satisfying 

(Ga,Gb)  =  8ab,  (r,gp)  =  s;,  <g‘fig  b)  =  8ab. 

Metric  tensors  are  introduced  on  each  configuration,  obeying  the  relations 

Gab  f  f  i  f "  B  ■  gxfS  it :  if .  ■  Sab  4.  it-  ■ 


(3) 

(4) 

(5) 

(6) 


The  notation  G  =  det  G,  g  =  det  g,  and  g  =  det  g  denotes  determinants  of  metric  tensors  given  in 
(6).  Since  we  often  use  volume  averaging  operations  over  the  crystal  element  to  define  certain 
tensorial  quantities,  we  restrict  the  basis  vectors  (3),  basis  covectors  (4),  and  metric  tensors  (6)  to 
be  constant,  but  not  necessarily  Cartesian,  within  each  element  (i.e.,  vie{,  Vmt,  and  rcur)  in  each 
configuration,  such  that  covariant  and  partial  differentiation  are  equivalent  operations.  However, 
these  variables  are  permitted  to  vary  from  volume  element  to  volume  element  (if  curvilinear 
coordinates  are  useful),  and  also  from  configuration  to  configuration. 

The  configurations  of  the  crystalline  volume  element  shown  in  Fig.  1  are  now  defined  in  turn. 
The  reference  configuration  volume  rref  consists  of  the  crystalline  lattice  as  it  existed  prior  to 
application  of  external  forces  (i.e.,  at  t  =  0),  such  that  it  is  free  of  traction  along  the  external 
boundary  sref.  It  may  or  may  not  contain  dislocations,  internal  residual  elastic  lattice  strains,  or 
residual  plastic  deformation.  The  current  configuration  volume  element,  rcur,  is  the  elastoplasti- 
cally-deformed  crystal  material,  with  possibly  non-vanishing  traction  vector  t  per  unit  reference 
area  applied  on  external  boundary  scur. 

The  local  deformation  gradient  f  for  material  points  with  local  coordinates  XA  within  the 
volume  element  is  defined  as  the  tangent  mapping 

f=T<Px  =  ^gaXGA.  (7) 


The  volume-averaged  deformation  gradient  F  for  the  crystal  element  is  then  defined  by  the  motion 
of  its  external  boundary,  which  is  equivalent  to  the  volume-averaged  local  deformation  gradient 
upon  invocation  of  Gauss’s  theorem  [61]: 


Fa 


1 

^ref 


X  (d£ref 


1 

^ref 


dxa 

OX* 


d^ref  — 

^ref 


(8) 


In  Eq.  (8),  (d.vrer)4  =  («rcf  )  ^ d.vrer  is  an  oriented  differential  surface  element  of  .vref  (Fig.  1),  with  nref 
a  unit  covariant  vector  normal  to  .vref. 
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Fig.  2.  Configurations  for  crystal  volume  element. 


An  average  elastic  stretch  tensor  Ve  is  associated  with  the  average  external  stress  applied  to  sw- 
The  intermediate  configuration  bml  (Fig.  2)  of  the  crystal  element  reached  upon  hypothetical 
instantaneous  elastic  unloading  from  the  current  configuration  via  the  inverse  of  the  elastic  stretch 
Ve~'  corresponds  to  null  traction  conditions  on  the  external  boundary  of  the  crystal  volume  el¬ 
ement  Dint  (i.e.,  the  traction  tint  =  0  along  as  shown  in  Fig.  2.  The  left  elastic  stretch  tensor  Ve 
is  determined  explicitly  from 


Ve  = 


(9) 


where  x  are  the  local  coordinates  of  the  external  boundary  of  the  element  corresponding  to 
traction-free  intermediate  configuration  bmt.  Configuration  bmt  arises  from  instantaneous  removal 
of  traction  along  the  boundary  of  Dint,  constrained  in  such  a  way  that  the  global  rotation  of  the 
volume  element,  Re  does  not  occur  upon  stress  relaxation.  Upon  unloading,  plastic  rear¬ 
rangements  are  idealized  as  rate  independent  and  inertial  effects  are  neglected. 

The  total  “elastic”  rotation  tensor  Re  is  determined  from  the  solution  of  the  following  integro- 
differential  equation  for  the  elastic  spin  and  associated  initial  conditions: 

RcRe  1  =—  [  fere-‘ dr ref,  Rc(?  =  0)  =  re(?  =  0)  =  1,  (10) 

Uef  J»ref 


where  1  is  the  identity  map  and  rc  is  the  local  elastic  and  rigid  body  rotation  exhibited  by  drref  as 
it  is  deformed  to  its  current  representation  drcur-  Note  that  RcT  =  Re_1  follows  from  (10)  and 
rcT  =  rc  1 ,  since  volume  averaging  preserves  the  anti-symmetric  property  of  the  spin. 

The  intermediate  configuration  buU  is  defined  by  the  net  unloading  procedure  Fe  1  =  ReTVe  ', 
i.e.  unloading  by  the  external  forces  and  subsequent  rotation  by  the  inverse  (i.e.,  transpose)  of  the 
average  lattice  rotation.  After  this  unloading  procedure,  local  coordinates  xa  describe  the  positions 
of  particles  within  the  volume  element  fint.  The  local  residual  deformation  gradient  f  is  the  tangent 
mapping 
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i=Tcpx  =  ^-A^®G\  (11) 

with  xf  =  (p'{x“{XA ,t),t)  the  motion  for  differential  volume  elements  between  bIC f  and  bml.  The 
volume-averaged  residual  deformation  gradient  F  for  the  crystal  element  is  deduced  from  the 
intermediate  configuration  coordinates  of  its  external  boundary,  and  is  equivalent  to  the  volume- 
averaged  residual  local  deformation  gradient  upon  invocation  of  Gauss’s  theorem  (compare  with 
(8))  “ 


_  IT  1  r  1  r  _ 

Fa  =  —  x«{dsrel)A=—  —  diyef  =  —  /  fA  dtVef-  (12) 

^ref  J sre f  ^  rel  J  sre{  Vre{  J  sre[ 

We  now  make  two  additional  assumptions  regarding  the  kinematics  of  single  crystalline  elas- 
toplasticity  on  a  “pointwise”  basis,  for  each  differential  volume  element  [62] 

f  =  fcfp,  f  =  fefP.  (13) 

In  Eq.  (13)i,  fc  =  vere  is  the  total  lattice  stretch  (ve)  and  rotation  (re)  for  a  differential  volume 

element  deformed  to  the  current  configuration,  and  fp  is  the  remaining  (plastic)  deformation 

gradient  attributed  to  the  history  of  motion  of  dislocations.  In  Eq.  (13)2,  fe  embodies  the  residual 
elastic  lattice  stretch  and  rotation  for  a  differential  volume  element  contained  within  the  externally 
unloaded  element  vmU  and  fp  is  the  residual  plastic  deformation  gradient  arising  from  dislocation 
motion.  In  order  to  illustrate  (13),  we  label  differential  volume  elements  as  drref  c  t'rer,  d5int  C  Smt, 
and  d/;cur  c  vcur  in  Fig.  3  with  differential  elements  visually  enlarged  relative  to  their  surrounding 
crystal  volume  elements  for  clarity.  The  local  deformation  gradients  of  (13)  then  act  as  mappings 
between  tangent  spaces  of  local  “differential  configurations”: 

f  :  T(drref)  -  T(dvCUI),  fp  :  r(diw) ->  T(drp),  fc  :  T(drp)  -  T(drcur), 

f  :  r(drw)  -  r(dgint),  fp  :  T(drref)  T(drp),  fc  :  T(di>p) T(dvmt), 


Fig.  3.  Tangent  mappings  and  volume  elements  at  multiple  length  scales. 
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Notice  that  of  the  six  tangent  mappings  in  (14),  only  f  and  f  are  necessarily  compatible  defor¬ 
mation  gradient  fields,  in  agreement  with  definitions  (7)  and  (11). 

The  total  plastic  deformation  gradient  FP  for  the  volume  element  emerges  from  the  following 
integro-differential  equation  written  in  terms  of  the  plastic  velocity  gradient: 

FPFP_1  =  —  [  fpfp-'  dr  ref,  (15) 

vie(  J  t>ref 

with  initial  conditions  FP(t  =  0)  =  fp(t  =  0)  =  1  if  a  perfect  reference  lattice  is  assumed.  Notice 
that  if  the  local  plastic  flow  is  isochoric,  i.e.  tr(fpfp  ')  =  0,  then  the  volume-averaged  plastic  flow 
of  (15)  is  volume-preserving  as  well,  i.e.  tr(FPFP  )  =  0. 

Next  we  assume  a  three-term  multiplicative  decomposition  for  the  total  deformation  F: 

f  =  fcfF=ff)  (16) 

T 

where  the  net  residual  deformation  gradient  F  enters  the  decomposition  as  shown,  in  accordance 
with  our  previous  definitions  of  Fe  and  F.  Combining  Eqs.  (8),  (12),  (13)2,  and  (16)  leads  to  the 
definition  of  F1: 

F‘  =  FFP~‘  =  —  (  [  fcfpdrref^)FP_1.  (17) 

Cef  \JvK f  / 

The  F1  (two-point)  tensor  can  be  thought  of  as  an  indicator  of  residual  elasticity  in  configuration 
bint  (and  corresponding  residual  stresses  and  interaction  energies).  The  deformation  tensor  F1 
contains  both  residual  elastic  (fe)  and  plastic  (fp,FP)  contributions.  It  is  clear  from  Eq.  (17)2  that 
F1  — »•  1  as  fc  — >  1  and  fp  — ■»  FP,  as  would  be  the  case  for  homogeneous  deformation  of  the  entire 
crystalline  element,  such  that  external  unloading  by  Ve  1  relieves  all  local  internal  stresses  and 
plastic  rearrangements  do  not  occur  upon  instantaneous  unloading.  On  the  other  hand,  if  the 
elastoplastic  deformation  fields  are  heterogeneous  throughout  the  volume  element,  F1  and  com¬ 
mensurate  residual  stresses  may  not  vanish.  Such  is  the  case  when  cellular  or  laminar  dislocation 
substructures  misoriented  from  one  other  by  finite  lattice  rotations  (embodied  here  locally  by  the 
orthogonal  part  re  of  fc)  evolve  and  refine  during  large  strain  deformations  in  FCC  metals  at  room 
temperatures  [43,44],  When  the  volume  element  encompasses  an  entire  single  crystal  (6ref  =  La  in 
Eq.  (1)),  we  thus  conclude  that  F1  represents  the  contribution  of  grain  subdivision  processes  to  the 
total  deformation  gradient  F  for  the  crystal,  as  originally  proposed  by  Butler  and  McDowell  [50], 
Additionally,  if  heterogeneity  of  local  deformation  and  stresses  due  to  elastoplastic  incompati¬ 
bility  (e.g.,  increased  slip  system  activity  and  stress  concentrations)  are  intensified  in  the  vicinity  of 
misoriented  grain  and  subgrain  boundaries  (cf^ experimental  and  numerical  results  in  [62]),  then 
one  would  expect  the  largest  contributions  to  F1  to  come  from  these  regions. 

The  deformation  tensor  F1  represents  in  an  average  sense  the  incompatibility  of  the  local  mi¬ 
croelastic  deformation  f  1  within  the  volume  element.  If  the  local  elastic  unloading  fe  1  is  uniform 
(and  hence,  compatible)  throughout  vcm,  then  fe  =  Fe,  fe  =  1,  and  F  =  f  =  FP  =  fp  =  fp,  such  that 
F1  =  1  (Fig.  3).  However,  F1  tells  us  nothing  about  the  compatibility,  or  lack  thereof,  of  the 
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average  recoverable  elastic  deformation  Fe  1  from  volume  element  to  volume  element.  If  we  re¬ 
gard  each  vre{  as  an  entire  single  crystal,  then  the  incompatibility  of  Fe  1  measures  the  inter¬ 
granular  incompatibility  between  grains,  while  if  we  regard  each  rref  as  a  subgrain,  then  the 
incompatibility  of  Fc  1  is  an  intragranular  measure.  By  incompatibility  of  Fe  1  we  mean  lack  of 
continuous,  single-valued  coordinates  x1  spanning  the  union  of  i  local  intermediate  configuration 
volume  elements  Vml  =  Urj^.  Thus,  while  the  x“  are  assumed  to  be  available  and  differentiable 
within  each  local  volume  element  v^v  they  may  not  be  so  in  the  global  configuration  Bmt  =  Vint.  If 
the  x“  are  multi-valued  or  discontinuous  in  the  global  configuration,  they  are  typically  called 
anholonomic  coordinates  [63]. 

We  now  appeal  to  the  continuum  theories  of  continuously  distributed  dislocations  [64-66] 
to  characterize  the  incompatibility  of  global  configuration  Bmt  (i.e.,  the  anholonomicity  or  lack 
of  integrability  of  Fe_1).  The  coefficients  of  the  “crystal  connection”  (cf.  [31,67])  acquire 
the  following  form  when  referred  to  the  current  configuration: 


r  =  Fe 

be  -a 


dFe 


0x* 


=  pe  pe-1  =  —pe-1  f 


a, 6’ 


(18) 


with  the  subscripted  comma  denoting  partial  differentiation  with  respect  to  spatial  coordinates  as 
indicated.  In  (18)  we  require  F'e  1  to  be  spatially  differentiable  to  first-order;  thus,  for  the  present 
discussion,  local  misorientations  across  grain  and  subgrain  boundaries  (i.e.,  between  volume 
elements  rCUr)  are  envisioned  as  steep  gradients  of  lattice  rotation,  as  opposed  to  actual  discon¬ 
tinuities  in  the  lattice  arrangement.  The  torsion  of  the  connection  r°bc  is  written 


t  =  -2F^F^ga 


g  ®g 


=  -2^7c]i 


(19) 


where  we  have  replaced  partial  differentiation  (subscripted  comma)  with  partial  covariant  dif¬ 
ferentiation  with  respect  to  the  Levi-Civita  connection  on  Bcm  whose  Christoffel  symbols  stem 
from  the  components  of  the  metric  gah  of  Eq.  (6)3  and  thus  are  symmetric  in  covariant  indices.  We 
use  the  notation  of  vertical  bars  for  such  covariant  differentiation.  Additionally,  the  bracketed 
indices  are  anti-symmetrized  according  to  2A[ab]  =  Aah  —  Aha. 

We  now  relate  the  torsion  tensor  t  of  Eq.  (19)  to  the  anholonomicity  of  the  inverse  elastic 
deformation  gradient  field  Fe  1 .  The  incompatibility  of  regions  within  global  configuration  BmX  is 
associated  with  the  net  Burgers  vector  B.  defined  in  terms  of  the  closure  failure  of  the  line  integral 
of  dx  =  Fe~'  dx  over  a  closed  loop  c  in  BmX  [65,66] 


B  =  — 


dx 


Feldx. 


(20) 


Applying  Stokes’  theorem  (cf.  [4]),  we  express  Eq.  (20)  as 


F^dx“=  I  F^dx“Adxh  = 


FhcFeJ  ncda, 


(21) 


where  dx"  A  dx*  =  sabcnc  d a  =  ffgeahcnc  da  is  the  differential  area  two-form  corresponding  to  ori¬ 
ented  area  element  ndu  bounded  by  curve  c  (Fig.  4),  which  in  turn  is  the  current  configuration 
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Fig.  4.  Burgers  circuit  in  current  configuration. 


image  of  Burgers  circuit  c.  The  wedge  product  is  denoted  by  A,  and  components  of  the  contra- 
variant  permutation  tensor  (i.e.,  alternator  tensor)  in  the  current  configuration  are  denoted  by 
tfhc  =  (  ffg)eabc ,  with  eabc  the  standard  permutation  symbols. 

The  incompatibility  is  expressed  in  terms  of  the  torsion  of  the  crystal  connection  by  combining 
(18)  (21)..  i.e. 


Bx  =  -^  J  Fe-vXh dx"  A  dxh .  (22) 

By  defining  the  operator 

curl(-)  =  (0(-)/0x)  :  £,  (23) 

we  may  write  for  (21)2 


B  = 


curl(Fe  1 )  •  n  da 


Ae  •  ndu, 


(24) 


where  Ae  =  curl(Fe~')  and  the  vector  ndu  =  {na da)ga.  Invoking  Nanson’s  formula  iid a  = 
.7*  'n  •  Fc  da  (cf.  [35]),  we  can  transform  (24)  to  an  area  integral  over  an  intermediate  configuration 
region  a 


B  =  /  n  JreFe_1AeT  d  a=  In-  Aeda, 


Ae 


(25) 


where  the  dot  product  in  the  integrand  of  (25)  denotes  contraction  of  the  covariant  index  of  n  with 
the  contravariant  index  of  Fe  '.  Note  that  Ac  only  measures  GNDs  and  does  not  account  for 
SSDs  (e.g.,  closed  dislocation  loops)  within  vie{  [68,69], 

Discussion  regarding  the  length  scales  inherent  in  the  line  integral  (20)  is  now  in  order.  Two 
distinct  length  scales  are  included  in  the  integral  (20)  and  are  intimately  related.  One  is  the  choice 
of  size  of  c:  we  let  D(c)  denote  the  diameter  of  a  circle  equivalent  to  c.  The  other  is  the  choice  of 
£ie f,  which  defines  the  size  of  the  volume  element  vre{  =  £]e{,  and  thus  rint,  and  implicitly  affects  the 
definition  of  the  average  elastic  deformation  gradient  for  the  element,  Fe.  Eq.  (20)  is  really  valid 
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only  for  D(c)  f  lre f,  if  we  take  fref  to  be  a  characteristic  diameter  of  the  volume  element  fmt 
considered.  Under  this  condition,  Fe  and  its  gradient  are  defined  in  regions  within  circuit  c.  If  the 
converse  is  true,  such  that  D(c )  <  £re{,  we  need  to  consider  the  incompatibility  of  the  inverse 
microelastic  deformation  within  the  volume  element,  fe  ',  instead  since  Fe  1  is  not  defined  within 
the  volume  element.  Repeating  the  derivation  leading  up  to  (24)  for  the  microelastic  deformation 
f6-1,  we  arrive  at  an  expression  for  a  local  Burgers  vector  be 

be  =  —  (j)  fe_1dx  =  f  curl(fc_1)  -ndu.  (26) 

Notice  that  bc  =  hc'gy  G  7’(drp)  resides  in  a  different  tangent  space  than  the  vector  B  of  (25).  We 
use  barred  Greek  indices  to  denote  components  in  this  space,  which  arises  from  the  local  plastic 
deformation  fp,  as  shown  in  Fig.  3.  Under  certain  conditions  we  can  further  characterize  (26)  in 
terms  of  the  spatial  gradient  of  the  inverse  of  the  residual  elastic  deformation  gradient  fc.  Since 
both  dt)int  and  d/;cur  inhabit  local  volume  elements  deformed  in  a  compatible  fashion,  we  can  define 
a  local  compatible  tangent  mapping  between  them:  f0-1  =  || :  r(drcur)  — >  T{ dSint).  If  no  additional 
local  plastic  deformation  occurs  upon  external  unloading  of  the  volume  element,  in  the  process 
symbolized  by  Fe  then  microscopically  we  have  fp  =  fp  and  drp  =  dvp,  as  shown  in  Fig.  5. 
Furthermore,  we  may  write 

f  =  fcfc,  (27) 

thereby  decomposing  the  total  microelastic  deformation  into  a  residual  part  fc  and  a  compatible 
part  fc  due  to  the  applied  stress  [70].  Substituting  (27)  into  (26)  and  invoking  the  compatibility 
conditions  /Pu  /[/i  =  0  then  gives 


Fig.  5.  Tangent  mappings  and  volume  elements  when  fp  =  fp. 
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demonstrating  that  bc  of  Eq.  (26),  representing  the  anholonomicity  of  the  local  elastic  deforma¬ 
tion  gradient,  depends  only  upon  the  spatial  gradient  of  the  (inverse  of  the)  residual  elastic  lattice 
deformation  map  fe.  and  not  higher  order  gradients  of  the  locally  compatible  and  recoverable 
elastic  lattice  deformation  field  fc. 


3.  Multiscale  balance  laws 


Neglecting  body  forces,  let  the  crystal  element  be  subjected  to  arbitrary  but  self-equilibrating 
surface  traction  t,  measured  per  unit  area  in  the  reference  configuration,  with  corresponding 
outward  unit  normal  covector  nref  (Fig.  2).  Let  sA“  represent  the  contravariant  components  of  the 
local  two-point  nominal  stress  tensor  (the  transpose  of  the  first  Piola-Kirchhoff  stress  tensor).  The 
traditional  (non-polar)  local  balances  of  linear  and  angular  momentum,  respectively,  are  written 
for  quasi-static  conditions  as 

<;  =  o,  faAsAb  =  sAafbA ,  (29) 

with  nrefAs  °  =  t“  on  Sref,  and  with  the  vertical  bar  denoting  covariant  differentiation  with  respect 
to  the  symmetric  Levi-Civita  connection  on  bK(.  Notice  that  Eq.  (29)  are  applicable  locally,  for 
points  (i.e.,  differential  elements)  within  the  crystal  element.  We  next  define  the  average  nominal 
stress  tensor  for  the  volume  element  as  (cf.  [71]) 

S  =  — —  f  sdrref  =  —  f  Xtgitdsvef,  (30) 

Cef  J  t>ref  Cef  J  sTcf 


where  Gauss’s  theorem  has  been  used  for  Eq.  (30)2,  along  with  the  assumptions  of  quasi-static 
conditions,  Eq.  (29)i,  and  stress  continuity  within  the  volume  element.  The  definition  of  the 
contravariant  Cauchy  stress  X  then  follows  as 


Iab  =  J-lFaASAb, 


(31) 


with  J  =  detF y/g/G.  The  macroscopic  rate  of  the  deformation  gradient,  F,  is  given  by 

i  =  —  [  f  d('rer  =  4  (  —  f  f  diyef  ]  =  —  [  x  (8)  dsVef, 

Cef  JvIer  dt\VIe{JViet  J  VIef  JSre[ 


(32) 


where  x(X,  t)  =  ^(x°gfl)  is  a  compatible  (material)  surface  velocity.  Consider  the  following 
equality  [72]: 


1 

Cef 


1 

^ref 


(x  -  FX)  (g)  (nref(s  -  S) )  dv 


ref- 


'  5ref 


(33) 


The  right-hand  side  of  Eq.  (33)  is  zero  for  linear  displacement,  constant  traction,  and  certain 
periodic  boundary  conditions  specified  with  respect  to  the  reference  configuration.  Assuming 
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henceforward  such  boundary  conditions  apply  as  representative  for  our  crystal  element,  the 
volume-averaged  stress  work  rate  per  unit  volume  is  found  as 


1 

^ref 


t  §a/?X  C^ref 


1 

^ref 


SAagabfhA  drref  =  SAaF 


■A) 


'  ^ref 


(34) 


where  we  have  used  Gauss’s  theorem  in  (34)  i. 

The  balance  of  energy  for  a  purely  mechanical  process,  considering  the  localized  form  of  (34) 
in  the  absence  of  heat  flux  due  to  conduction  or  internal  heat  supply,  is  written  as 

uo  =  SAjaA,  (35) 

with  m0  the  local  internal  energy  per  unit  volume  in  the  reference  configuration.  The  corresponding 
averaged  energy  balance  then  follows  from  Eq.  (34)2  and  (35)  as 

Uo  =  SAaFaA,  (36) 

where  U0  =  (rrer)  1  /'  m0 duref  is  the  average  internal  energy  per  unit  reference  configuration 
volume.  Notice  that  Eqs.  (35)  and  (36)  are  constitutive  assumptions,  prohibiting  conjugate  mi¬ 
crostresses  to  second-order  gradients  of  deformation  from  performing  mechanical  work  on  the 
external  boundaries  of  diyef  and  rref,  respectively.  These  rather  strong  assumptions  are  in  agree¬ 
ment  with  previous  treatments  by  Teodosiu  [3],  Steinmann  [32],  Acharya  and  Bassani  [33,73], 
Menzel  and  Steinmann  [9],  Bammann  [34],  Bassani  [74],  and  Regueiro  et  al.  [39],  On  the  other 
hand,  they  are  contradictory  to  the  models  of  Teodosiu  [1],  Dillon  and  Kratochvil  [5],  Naghdi 
and  Srinivasa  [38,28],  Le  and  Stumpf  [30,31],  Shizawa  and  Zbib  [16,17],  and  Gurtin  [13,36], 
Introducing  )/0  as  the  average  specific  entropy  per  unit  reference  volume,  the  Clausius-Duhem 
inequality  is  written 

no  >  0,  (37) 

again  in  the  absence  of  heat  flux  due  to  conduction  or  internal  heat  supply.  The  average  Helm¬ 
holtz  free  energy  per  unit  volume  in  the  reference  configuration,  i//0,  is  then  defined  by 

<Ao  =  U0-  rj06,  (38) 

with  6  the  absolute  temperature.  Substituting  Eqs.  (36)  and  (38)  into  the  Clausius-Duhem  in¬ 
equality  (37)  and  assuming  stationary  temperature  then  yields  a  reduced  form  for  the  entropy 
inequality: 

SAaF°A  >  K  (39) 

The  remaining  balance  laws  at  the  level  of  the  volume  element  are  now  presented: 
p0  =  pj,  S$  =  0,  FaASAb  =SAaFbA. 


(40) 
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Eq.  (40)i  provides  the  usual  relationship  between  current,  p,  and  reference,  p0,  average  mass 
densities.  Eq.  (40)2  and  (40)3  follow  from  volume  averaging  the  microscopic  balance  laws  in  Eq. 
(29)  and  by  assuming  that  boundary  conditions  imposed  on  the  crystal  volume  element  are  such 
that  the  right-hand  side  of  Eq.  (33)  is  zero  throughout  the  deformation  process.  The  last  of  (40) 
leads  to  symmetry  of  the  Cauchy  stress  E,  as  defined  in  Eq.  (31). 

Notice  that  inequality  (39)  is  written  in  terms  of  two-point  tensors  S  and  F  and  the  free  energy 
per  unit  reference  volume  i//0.  We  shall  later  find  it  useful  to  express  (39)  in  terms  of  quantities  with 
all  components  referred  to  intermediate  configuration  bim  (i.e.,  the  relaxed  intermediate  config¬ 
uration  of  elastoplasticity),  which  we  view  as  the  most  convenient  configuration  for  deducing 
thermodynamic  restrictions  and  posing  constitutive  assumptions  (see  also  [16,17,39]).  Expression 
of  (39)  in  the  intermediate  configuration  is  achieved  here  by  first  mapping  to  the  current 
configuration,  then  pulling  back  by  the  two-point  elastic  deformation  gradient  Fe.  From  (31)  we 
have 


SAaFaA  =  JY?h  gbcFcAFflA  =  .FL“hL 


bai 


—Lba 


(41) 


with  L  the  spatial  velocity  gradient.  Substituting  into  (39)  and  multiplying  by  the  inverse  of  the 
Jacobian  invariant  J  =  det  F \Jg/G  >  0  then  gives  the  Clausius-Duhem  inequality  in  the  current 
configuration,  i.e. 


(rhLha  +  ilr0j-l)>iii,  (42) 

where  the  free  energy  per  unit  current  volume  \sij/  =  J  From  the  symmetry  of  E  and  g  ',  and 
using  the  standard  identities  J  lJ  =  0  and  j  =  JgahLha,  Eq.  (42)  can  be  simplified  to  read 


(43) 


with  f/\  the  Lie  derivative  with  respect  to  the  spatial  velocity  field  v  =  xof',  and  with 
(A,  B)  =  tr(AB)  the  dual  product  defined  for  second-rank  tensors  A  and  B.  We  next  define  the 
velocity  gradient  L  referred  to  intermediate  configuration  bm[  as 


L%  =  Fe-r 


Ta  Fe 


(44) 


and  we  define  the  mixed-variant  elastic  second  Piola-Kirchhoff  stress  S  and  the  Mandel  stress  [75] 
M  as 


S*f>  =  J*FZ 


^gcbFl 


M°p  =  jsf: 


z acgcbF\ \ 


—  re  v* 

- 


(45) 


where  Ce“  =ga^Fe"gabF^  and/6  =  det F c\/g/g.  Then  from  relations  (39),  (41)  and  (45)2,  we  can 
obtain 


(46) 
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Multiplying  (46)  by  J  1  >  0  and  using  the  relation  J-1  =  —J  2J  =  —  J~l F*aF~iA ,  along  with 
i jt  =J  V0,  then  gives 


-  ij/F* 


which  is  the  reduced  entropy  inequality  mapped  to  configuration  bmi,  with  energetic  quantities 
defined  on  a  per-unit-volume  basis  with  respect  to  local  element  volume  fmt. 


4.  Crystal  plasticity  model 

Hereafter  we  more  fully  develop  a  continuum  crystal  plasticity  model,  focusing  on  a  formu¬ 
lation  couched  at  a  single  length  scale  £ief  rather  than  one  phrased  in  terms  of  volume  averages. 
This  is  necessary  for  typical  numerical  implementations,  wherein  each  volume  element  /.'rei  cor¬ 
responds  to  the  local  region  about  an  integration  point  in  a  finite  element  simulation,  for  example. 
The  microscopic  deformation  gradients  f,  f,  fe,  fp,  fc,  fp,  and  fc  were  introduced  in  Section  2  for 
illustrative  purposes  and  are  not  calculated  explicitly  in  this  model.  The  same  can  be  said  for  the 
microscopic  kinetic  variables — such  as  t,  s  and  u0,  for  example — of  Section  3.  Instead,  we  consider 
her  (y  the  evolution  of  the  average  deformation  gradients  for  the  crystalline  volume  element  (e.g.,  F, 
Fe,  F1,  and  FP),  the  average  stresses  (e.g.,  S  and  2),  and  the  average  energies  (e.g.,  U0  and  i j/0).  Since 
the  microscopic  deformations  and  local  lattice  rearrangements  within  the  volume  element  are  no 
longer  explicitly  tracked,  we  are  unable  to  invoke  Eq.  (9)  to  determine  Ve,  Eq.  (10)  for  Re,  Eq.  (15) 
for  FP,  Eq.  (17)  for  F1,  or  Eq.  (26)  to  calculate  the  local  incompatibility  bc  within  microscopic 
subregions  of  the  crystal  element.  Likewise,  we  are  unable  to  calculate  local  residual  stress  fields 
and  elastic  energies  associated  with  fe  in  the  model  forthcoming  in  Section  4.  Instead,  we  must  rely 
upon  additional  constitutive  assumptions  and  corresponding  balance  relations  to  ready  our  model 
for  implementation — these  are  considered  in  detail  in  what  follows.  Of  course,  the  constitutive 
assumptions  we  make  hereafter  are  motivated  by  the  physically-based,  multiscale  averaging 
treatment  of  Sections  2  and  3. 

4.1.  Kinematics,  balance  laws,  and  thermodynamics:  summary 

The  fundamental  thermomechanical  relations  already  derived  in  Sections  2  and  3  and  applied 
now  to  a  “single-scale”  crystal  plasticity  formulation  are  restated  here  for  clarity  and  ease  of 
reference. 

Deformation  gradient: 

„  _  _ _ 

F  =  ^  =  V^VR^F".  (48) 


Dislocation  density  tensor  (GNDs): 
Ae=.7cFe  1  (curlFc 


(49) 
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Mechanical  stresses: 


M  =  g-'CcS  =  JeFeT2gFe_T  =  J_1FeTFSgFe~T. 
Mass  conservation: 


(50) 


Po  =  pJ- 


(51) 


Balance  of  linear  momentum  (quasi-static): 

Div(S)  =  0.  (52) 

Balance  of  angular  momentum  (quasi-static,  intermediate  configuration): 

S  -  ST  =  0.  (53) 

Balance  of  energy  (mechanical  case,  intermediate  configuration): 

(MT,r_]FF_^T)  =  J-'U0.  (54) 

L 

Entropy  inequality  (mechanical  case,  intermediate  configuration): 

(55) 


4.2.  Free  energy  and  consequences  of  the  dissipation  inequality 

We  make  the  following  general  assumption  regarding  the  dependency  of  the  Helmholtz  free 
energy  function  for  the  crystal  volume  element,  referred  to  relaxed  intermediate  configuration  bmt 
and  neglecting  temperature  effects: 


f  =  f(Ce,\',Ae,sss,g) 


-  P^SSig«p)- 


(56) 


The  covariant  elastic  strain  tensor  C! g  =  FfgahFf  is  included  to  model  the  change  in  average 
energy  with  a  change  in  external  loads,  a  standard  assumption  in  finite  crystalline  elastoplasticity 
theories  (cf.  [76]).  The  left  stretch  tensor  derived  from  F1,  denoted  by  V1,  is  included  to  reflect  the 
contribution  to  the  free  energy  from  residual  microelasticity  within  the  volume  element,  and  is 
non-negligible  when  the  deformation  within  the  volume  element  is  heterogeneous  (e.g.,  during 
grain  subdivision).  Kratochvil  [77]  and  Lion  [78]  made  similar  constitutive  assumptions.  The 
elastic  energy  due  to  the  average  elastic  curvature  in  the  volume  element  attributed  to  GNDs  is 
reflected  by  the  inclusion  of  Ae,  a  particular  assumption  also  suggested  explicitly  by  Bammann 
[34],  Gurtin  [36],  and  Regueiro  et  al.  [39],  Since  the  GND  density  tensor  does  not  include  a 
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measure  of  the  total  length  of  all  dislocation  lines  (e.g.,  SSDs  consisting  of  closed  dislocation 
loops  and  dislocation  lines  of  opposing  signs),  the  scalar  parameter  ess  =  b^/J> ^  is  included  to 
model  the  net  contribution  of  the  elastic  self-energy  of  the  SSDs  to  the  total  free  energy  (cf. 
[79,80]).  The  scalar  norm  of  the  Burgers  vector  (e.g.,  the  lattice  parameter)  is  denoted  by  b ,  while 
the  line  length  of  SSDs  per  unit  intermediate  configuration  volume  is  written  pss.  We  regard  eSs  as 
a  local  residual  lattice  strain  measure  due  to  the  presence  of  SSDs  [34]. 

It  is  emphasized  that  all  independent  variables  in  the  free  energy  function  are  “elastic”  vari¬ 
ables,  in  contradiction  to  works  of  Naghdi  and  Srinivasa  [38,28],  Gurtin  [13],  and  Svendsen  [37], 
among  others,  who  include  a  dependency  upon  the  plastic  deformation  gradient  and/or  higher 
gradients  of  plastic  deformation.  We  treat  the  plastic  deformation  gradient  Fp  as  a  continuum 
idealization  of  the  rigid  sliding  of  portions  of  the  lattice  due  to  relative  motion  of  elastic  blocks  of 
material  contained  between  slip  planes.  Such  motions  clearly  do  not  alter  the  energetic  properties 
of  material  within  the  blocks,  and  therefore  should  not  influence  the  Helmholtz  free  energy. 
Reasoning  behind  inclusion  of  the  metric  g  will  be  explained  later. 

Notice  that  each  of  the  constitutive  variables  in  (56)  is  invariant  with  respect  to  superposed 
rigid  body  motions  in  the  current  configuration:  x  — >  Qx  +  c,  with  Q  €  S03  a  rigid  body  rotation 
matrix  and  c  a  constant  translation  vector: 


Fc  -»  QFe  =>■  Ce  -►  Ce,  F‘  F‘  V1  — >  V1, 

Ae  =/eFe-1(curlFe^1)T  — ^/e(detQ)Fe-1QT(curl(Fe-1QT))T  =  Ae,  (57) 

Pss  Pss  =>■  £ss  — >  ess,  g  ~ >  g- 

Even  though  F1  is  unaffected  by  rigid  Lra n s I'o rm at i o n s  o f the  current  configuration  [77],  we  choose 
not  to  include  the  rotational  part  of  F1,  denoted  by  R1,  in  the  free  energy  function  because  rigid- 
body  rotations  of  the  volume  element  in  the  unloaded  configuration  bint  do  not  influence  the 
stored  elastic  energy. 

We  define  the  conjugate  stresses  to  the  independent  state  variables  included  in  (56): 


Se 


~i  __  Si ji  ~GN  _  8i ji 

~  SV  ’  “  8Ae  ’ 


sss=d^ 


8& 


ss 


(58) 


Notice  that  SGN  is  a  couple  stress,  with  units  of  force!  length,  while  the  remaining  stress  measures, 
Sc,  S1,  Sss,  and  Sg  possess  the  standard  stress  dimensions  of  for  cel length2 .  Expanding  f  in  (55) 
using  the  definitions  of  (58)  gives 


L^M 


—  X!/F  AFfl  >S 


c^  +  s\ 


?,  +  S™Ae 


+  SSS&SS  +  SS  gap. 


(59) 


where  the  dual  map  [81,82]  is  denoted  by  ()*  and  satisfies  (A“hf  A'h“.  From  the  three-term 
multiplicative  decomposition  of  the  deformation  gradient,  F  =  FeF‘FP,  of  Eqs.  (16)  and  (48),  the 
velocity  gradient  in  the  intermediate  configuration,  L,  can  be  partitioned  as 
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L  =  Fe-1  FF~^Fe  =  F-'F  +  FF1"1,  +  FFV  V"1, .  (60) 

L  Le  V  LP 

Further  algebraic  manipulations  give 

(L,Mt)  =  ^Sg-1,^Ce^  +  (MT,Li  +  Lp),  FFj  =  vr1  +  ViRR^VM,  (61) 

L‘  W* 

where  the  spin  W1  is  skew  via  g^W'y  =  —gxp  wf .  Using  Eqs.  (60)  and  (61)  in  inequality  (59)  then 
results  in 


Q  5“rltf-^)c;+((Mj- 


-  Ss%s  +  (M.f  -  >  0. 


(62) 


Assuming  that  (62)  must  hold  for  independent  specification  of  the  elastic  strain  rate  and  the  other 
rate  variables  [83,84]  we  arrive  at  standard  relationships  between  the  (mechanical)  elastic  stresses 
S,  E,  S,  and  T,  and  the  lattice  stress  Sc  conjugate  to  Cc 


2Se  =  2-^1-  =  Sir1  =  r F^EF6-*  =  J-'FSFe~*  =  J“1Fe~1TF*,  (63) 

0Ce 

where  T  =  ST  is  the  contravariant  first  Piola-Kirchhoff  stress  tensor,  and  the  final  equality  follows 
from  (40)3.  More  familiar  and  compact  constitutive  equations  for  the  Cauchy  stress  and  first 
Piola-Kirchhoff  stress  are  then  derived  readily  with  the  chain  rule  of  differentiation 


Zah  =  2 Je-1  Fe“Fe*  =2Je-1^-^  2  ^ 


QC 


dg 


ab 


dgab  ’ 


=SC^/6goi 


TaA  =  2 JF7lAF? 


ffp-  =  J  F-'A  a.  r  =  J  2L  „»  ~  '2k 

act,  "  ■■  err  err  err 


=0i A/0U 


(64) 


(65) 


The  final  equality  in  (64)  is  rigorous  only  when  o/e  '/eg  =  0  (e.g.,  Cartesian  current  coordinates), 
while  the  final  equality  in  (65)  is  ensured  only  when  0//0F  =  0  (e.g.,  when  the  total  inelastic 
deformation  is  isochoric).  We  next  define  (cf.  [85]) 
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P  =  Mt  -  \ji\  =  JeFe~  (Eg  -  i/d)  Fe  =  F  (FtTG  -  iA0l0)  F~!  (66) 

S - V - '  S - v - ' 

p  p 

as  the  push-forward  of  the  mixed-variant  reference  configuration  Eshelby  1  [86]  energy-mo¬ 
mentum  tensor  P,  or  equivalently  as  the  pull-back  of  the  current  configuration  Eshelby  stress 
tensor  p.  In  Eq.  (66),  10,  1,  and  1  denote  mixed-variant  identity  maps  on  configurations  bre{,  bmt, 
and  bcm,  respectively.  Upon  substituting  (63)  and  (66)  into  (62)  and  assuming  that  g  =  0,  we  arrive 
at  a  reduced  form  of  the  dissipation  inequality: 

(P V^1  -  S1* ,  V1  +  (V'pr1 ,  W1)  -  (SGNT,  Ae)  -  (Sss,  ess)  +  (P,  Dp  +  Wp)  ^  0.  (67) 

Notice  that_  both  the  plastic  deformation  rate  Dp^=g^(gLp),^  and  the  plastic  spin 
Wj  =  g/‘,(gLp)|(j/;  contribute  to  the  dissipation,  as  does  the  spin  associated  with  residual  elas¬ 
ticity/subdivision  W1.  Note  also  the  prominent  role  of  the  Eshelby  stress  tensor  P  as  a  force 
conjugate  in  the  plastic  dissipation  term  (see  also  [85,87]). 

4.3.  Constitutive  model 

A  more  explicit  form  of  the  strictly  mechanical,  intermediate  configuration  free  energy  function 
(56)  is  now  proposed  on  physical  grounds: 

^  =\ <Ee,  (Ceff(E‘))  :  Ee>  +  clfl{  E\  Ee)  +  c2n\J  (E\  E‘)  +  c3/r(eSs)2  +  PasMK,  A* ).  (68) 

The  strains  Ec  and  E1  follow  from  the  decomposition  of  the  total  strain  tensor  E  referred  to  bml: 

E  =  ^ F~T(FTF  -  10)F_1  =  Ec  +  E‘  +  Ep,  (69) 

where  the  elastic  (Ee),  heterogeneity  (E1),  and  plastic  (Ep)  parts  are  given  by 

Ec  =  ^(g_1Cc  —  i),  Ei  =  ^(l-yi-Tyi-1),  Ep  =  ^Fi_T(l  —  Fp_TFp_1)Fi_1.  (70) 

Notice  that  Ec  and  E1  are  functions  of  Ce  and  V1,  respectively,  in  agreement  with  the  general  form 
for  the  free  energy  in  proposition  (56).  The  mixed-variant  GND  variable  A°#  is  defined  as 

A#  =  Acg.  (71) 


1  Eshelby  [86]  in  fact  used  the  negative  of  our  P  as  the  referential  energy-momentum  tensor,  with  the  material 
gradient  of  displacement  in  place  of  the  deformation  gradient.  We  follow  here  for  convenience  the  definition  and  sign 
convention  of  Le  and  Stumpf  [87]  and  Le  et  al.  [76], 
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Components  of  the  fourth-rank  effective  elastic  modulus  tensor  for  the  volume  element  in  the 
intermediate  configuration  are  denoted  by  (Ceff)7„|,  and  p  is  a  characteristic  elastic  constant  for  a 
perfect  reference  lattice  (e.g.,  an  average  reference  shear  modulus).  Finally,  ci,c2,  and  c3  are  di¬ 
mensionless  scalar  constants,  and  /iGN  is  a  (possibly)  evolving  scalar  parameter  that  will  be  ad¬ 
dressed  more  later  (Eqs.  (72)  and  (73)),  with  dimensions  of  length2.  Notice  the  role  of  g  in  defining 
the  mixed-variant  tensors  Ec  and  Ae#,  rendering  its  inclusion  in  the  general  free  energy  function 
(56)  a  necessity  [87], 

We  now  discuss  the  physical  reasoning  behind  our  choice  of  each  term  in  the  specific  free  energy 
function  (68).  The  first  term,  \  (Ee,  (Ccff(E'))  :  Ee),  is  reminiscent  of  the  quadratic  form  seen  in 
finite  linear  hyperelasticity,  differing  only  in  the  sense  that  the  effective  elastic  modulus  tensor  for 
the  volume  element,  referred  here  to  the  intermediate  configuration  bmU  is  assumed  to  depend 
explicitly  upon  E1.  This  assumption  is  intended  to  reflect  the  influence  of  local  microelastic  ro¬ 
tations  which  arise,  quite  possibly,  from  grain  subdivision  and  intragranular  substructure  de¬ 
velopment  during  plastic  deformation  (cf.  Flughes  et  al.  [43],  Butler  and  McDowell  [50]).  For 
example,  for  a  single  crystal  consisting  of  subgrains  exhibiting  a  random  distribution  of  local 
misorientations  unrestricted  in  magnitude,  the  effective  moduli  will  approach  those  of  a  random 
polycrystal  (e.g.,  elastic  isotropy  in  many  metals). 

The  second  term,  cjilE'.  Ee),  accounts  for  the  aforementioned  amplification  of  internal  re¬ 
sidual  microstress  fields  (and  the  corresponding  internal  elastic  energy)  at  flexing  subgrain 

boundaries  [59,60]  with  increases  in  the  applied  stress.  The  third  term,  c2jJ.\J (E\  E1),  represents  a 

portion  of  the  stored  energy  of  cold  work  attributed  to  heterogeneous  elastoplasticity  within  the 
volume^  element.  The  linear  dependency  of  energy  upon  the  effective  “meso-incompatibility” 
strain  E1  is  motivated  from  previous  solutions  obtained  from  computational  micromechanics  [88] 
in  which  low-angle  and  high-angle  grain  boundary  arrangements  were  assigned  as  initial  condi¬ 
tions  to  a  deforming  polycrystalline  aggregate.  Since  the  present  model  is  intended  for  single 
crystals,  we  presume  that  data  for  low-angle  grain  boundary  arrangements  from  this  previous 
study  are  most  applicable  (Fig.  6(a)),  and  our  form  of  the  free  energy  term  reflects  an  assumption 
of  similitude  between  single  crystals  and  polycrystals  with  low-angle  boundaries.  Thus,  this  term 
should  be  regarded  as  a  rough  initial  approximation. 

The  fourth  term,  /?GNp( Ae#,Ae#),  represents  a  quadratic  dependency  of  the  free  energy  upon  the 
average  density  tensor  of  GNDs.  Such  an  assumption  is  rather  standard  in  gradient  finite  strain 
single-  and  polycrystalline  elastoplasticity  theories  in  the  literature  (cf.  [9,32,34,37,39]),  as  is  the 
linear  dependence  upon  the  total  length  of  SSDs  per  unit  volume,  pss  =  (sss/b)2  [34,39,37].  Ad¬ 
ditional  experiments  and/or  numerical  simulations  are  clearly  needed  to  determine  the  parameters 
ci,  c2,  and  c3  for  a  particular  material. 

We  remark  that  the  parameter  /?GN  provides  a  squared  effective  “radius  of  non-local  action”  in 
our  model  (see  also  later  Eqs.  (90)  and  (93)).  As  /?GN  — >  0,  the  contribution  of  GNDs  to  the  free 
energy  function  becomes  negligible,  and  the  model  assumes  a  local  character.  We  hypothesize  that 
fref  =  v^ef — the  size  of  the  representative  crystal  volume  element  (Eq.  (1)) — should  influence  the 
value  of  /?GN  in  some  way.  In  other  words,  the  size  of  the  local  crystal  volume  element  over  which 
the  average  dislocation  density  tensor  Ae  is  calculated  should  have  some  effect  upon  the  range  of 
non-local  interactions.  For  example,  if  we  assume  that  the  strength  of  non-locality  is  inversely 
related  to  scale  (cf.  [24]),  then  an  obvious  choice  is 
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(a)  l|(2/3)<Ei,E‘> 


(b)  Applied  effective  total  strain 

Fig.  6.  Residual  elastic  energy  versus  \J (2/3) (E1 ,  E1)  (a)  and  evolution  of  \J (E1,  E‘)/(Ee,  Ee)  versus  applied  effective 
strain  (b). 


P  GN  —  ( 


2 

GN 


(72) 


with  /:"GN  a  length  parameter  that  is  characteristic  for  a  given  material  (cf.  [39])  and  with  c4  another 
dimensionless  material  constant.  In  other  words,  if  we  regard  £GN  and  c4  as  fixed  parameters 
(constants)  for  a  given  material,  then  the  non-local  radius  will  scale  with  the  size  of  the  volume 
element  according  to  iK f4/2.  Of  course,  we  would  now  need  a  more  extensive  series  of  tests  to 
determine  the  two  values  £GN  and  c4,  rather  than  just  /1GN.  Interestingly,  if  we  take  fref  =  La,  the 
grain  size,  and  c4  =  1,  then  a  Hall  Petch  type  relation  [18,19]  for  hardening  due  to  GNDs  is 
acquired  via  Eq.  (72). 

Additionally,  one  may  argue  that  since  the  radius  of  non-locality  should  reflect  a  characteristic 
dimension  of  the  microstructure,  such  as  cellular  structure  size  in  ductile  single  crystals  that 
subdivide  (cf.  [7]),  then  the  effective  length  fGN  should  be  an  evolving,  rather  than  stationary, 
parameter.  In  our  framework,  we  could  easily  extend  (72)  to 
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Pgn  = 


^GN  —  •^Gn(V1), 


(73) 


allowing  the  effective  scaling  factor  fCN  to  evolve  with  the  heterogeneity  parameter  (i.e.,  the  grain 
subdivision  metric)  V1.  If  Eq.  (73)  is  invoked,  then  the  presence  of  f G n  in  the  free  energy  function 
(68)  will  influence  the  values  of  thermodynamic  force  conjugates  calculated  via  Eqs.  (58)  con¬ 
tributing  to  the  dissipation  inequality  (67). 

Consider  a  boundary  value  problem  where  F  and  its  rate  are  prescribed  incrementally  as  a 
function  of  time.  In  order  to  characterize  completely  the  state  of  the  material  from  the  kinematic 
standpoint  of  Eq.  (48),  we  need  to  specify  the  deformation  gradient  measures  Fc,  F1,  and  FP.  In  a 
hyperelastic  setting  (cf.  [89,90]),  a  standard  approach  is  to  formulate  evolution  laws  for  (time  rates 
of)  FP  and  F1,  thus  leaving  Fe  to  be  determined  from  the  product  F(F‘FP)_1,  assuming  that  the  rate 
of  F  and  the  driving  forces  for  each  inelastic  deformation  measure  are  known  at  the  beginning 
of  each  time  increment.  Such  is  our  approach:  the  requisite  evolution  equations  for  FP  and  F1  are 
considered  next. 

The  time  rate  of  average  plastic  deformation  gradient  in  the  single  crystal  volume  element  is 
assumed  to  follow  from  the  standard  kinematic  framework  of  finite  crystal  plasticity  [80,91-93], 
i.e., 


4^ 

i=i 


y's'  <g)  m'  F, 


l p 


(74) 


where  y‘,  s',  and  m'  are  the  shearing  rate,  slip  direction  (a  unit  vector),  and  slip  plane  normal  (a 
unit  covector)  for  slip  system  i,  all  defined  with  respect  configuration  bp  of  Fig.  2.  As  written  in 
(74),  Lp  =  FPFP  is  the  mixed-variant  average  plastic  velocity  gradient  in  configuration  bp.  When 
the  single  crystalline  lattice  is  initially  heterogeneous  in  the  reference  state,  or  when  the  lattice 
deforms  heterogeneously  within  the  volume  element,  f,  s',  and  m'  are  understood  to  be  suitably- 
defined  spatial  averages  of  their  fluctuating  local  counterparts.  In  accordance  with  classical  crystal 
plasticity  theory,  the  representation  of  the  slip  directions  and  slip  plane  normals  is  assumed 
identical  in  each  of  the  unloaded  configurations  bK{  and  bp.  The  slip  direction  unit  vectors  and  slip 
plane  normal  covectors  are  orthogonal  in  the  unloaded  configurations  (i.e.,  (m',s')  =  0),  and  are 
typically  given  as  initial  conditions  in  a  boundary  value  problem.  However,  their  representation 
changes  in  the  current  configuration  as  a  result  of  the  elastic  deformation  gradient  Fe  (cf.  [76,93]) 


(Sr=Fe;5i(sf, 


(75) 


In  Eq.  (75),  s'  and  m'  are  push-forwards  of  their  counterparts  in  bp  (we  have  used  barred  Greek 
indices  for  components  in  bp)  and  are  not  necessarily  of  unit  length  when  the  elastic  stretch  Ve  is 
significant,  although  each  pair  of  s'  and  m'  does  remain  orthogonal  in  configuration  Z>cur,  as  is 
easily  verified  by  direct  calculation.  These  orthogonality  relations  effectively  prohibit  dislocation 
climb,  a  non-isochoric  process.  Two  useful  consequences  arise  from  kinematic  assumption  (74) 
and  the  slip  plane  normal-slip  direction  orthogonality  relations.  The  first  is 
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jp  =  yptr(Lp)  =  y'W,  s')  =  0,  (76) 

i 

meaning  that  the  plastic  flow  is  isochoric.  The  second  is  a  reduced  form  of  the  “purely  plastic” 
dissipation  in  the  entropy  inequality  (67): 

(P,  Dp  +  Wp>  =  (j^PF^LP)  =  (P,Lp)  =  (77) 

P  Z 

where  P  is  the  Eshelby  stress  tensor  mapped  to  “plastic”  configuration  bv  and 

(78) 

is  the  projected  Eshelby  stress_on  system  i.  We  could  equivalently  use  the  pull-back  of  the 
transposed  Mandel  stress,  F1_1MTF1_r,  in  place  of  P  in  (77)  and  (78)  since  the  plastic  flow  is 
volume-preserving.  By  assuming  a  particular  flow  potential  0  for  the  slip  rates  yl,  we  can  ensure 
that  the  dissipation  in  (77)  is  non-negative,  such  that  the  entropy  inequality  is  satisfied  auto¬ 
matically  when  the  contributions  of  F1,  Ae,  and  pss  to  the  dissipation  are  neglected.  We  conve¬ 
niently  prescribe 


0=^2&(p,f  \sgn,sss) = 

i=l  i 


K‘ 

m  +  1 


P_ 

k‘ 


m+ 1 


(79) 


where  each  slip  potential  is  a  scalar  function  of  its  arguments.  In  Eq.  (79)2,  k  >  0  (dimensions 
of  1  / time)  and  m  (dimensionless)  are  material  parameters  which  we  assume  are  constant  on  each 
slip  system,  while  k‘(F'.  SGN,  Sss)  >  0  are  scalar  slip  resistances  (dimensions  of  for  cel  length2)  that 
evolve  with  inelastic  deformation.  The  slip  rates  are  then  assumed  to  adhere  to 


0<2>‘ 

cp1 


ksgn(p‘) 


(no  sum  on  /), 


(80) 


a  power-law  form  reminiscent  of  typical  flow  rules  proposed  for  classical  crystal  viscoplasticity 
[94],  except  that  the  resolved  Eshelby  stress  is  used  as  a  conjugate  variable  to  the  slip  rates  in  our 
theory.  Substituting  (80)  into  (77),  we  find  that  the  purely  plastic  dissipation  is  unconditionally 
non-negative: 


X'Ei'wV) 


>  o. 


(81) 


Our  thermodynamic  motivation  for  choosing  the  resolved  Eshelby  stress  as  a  driving  force  for 
viscoplasticity  is  clear  from  (67)  and  (81).  Our  physical  reasoning  for  selecting  the  Eshelby  stress 
stems  its  relationship  to  the  configurational  force  on  a  dislocation  [85],  as  opposed  to  the  usual 
resolved  Cauchy  stress  of  Schmid’s  law  [95],  We  also  cite  the  work  of  Le  et  al.  [76]  who  used 
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thermodynamic  and  physical  arguments  similar  to  ours  to  motivate  usage  of  the  Eshelby  stress  as 
a  driving  force  for  plastic  deformation  in  a  finite  strain,  rate-independent  single  crystal  plasticity 
model. 

We  propose  the  following  dependency  for  the  positive  scalar  slip  resistances  k'\ 


k'(F\S  gn,Sss)  =  F'*(SGN)TFi,^£^)  =  k'(Sgn,Sss),  (82) 

gGN  5SS 

with  SG^  and  Sss  complete  pull-backs  of  SGN  and  Sss  to  configuration  bp,  and  with 
J1  =  det  F'^/det  g/detg  >  0  (g  is  now  introduced  as  a  local  metric  on  the  anholonomic  space  bp). 
The  thermodynamic  forces  SGN  and  S'88  satisfy 


/  oGN\  _  — 1  rpi*-ct  ^SS 

with 


0£ss  dess  ’ 


(83) 


Ae  =  J'F^A^F1-*,  ess  =  V7%s  (84) 

the  densities  of  GNDs  and  SSDs,  respectively,  pulled  back  from  Z?int  to  bp.  As  a  more  particular 
constitutive  assumption,  we  represent  the  hardening  variable  for  each  slip  system  as  a  linear 
combination  of  the  thermodynamic  force  conjugates  to  the  densities  of  GNDs  and  SSDs  parti¬ 
tioned  to  all  systems,  i.e., 

k  =k0+  \J  /I'gn^gn/’gn  +  b'SSjSJss,  (85) 

where  k'0  is  the  initial  friction  or  threshold  stress  for  system  i,  tiGN/.  and  h'SSj  comprise  n  x  n  slip 
system  interaction  matrices  (dimensionless  units)  to  be  obtained  from  experimental  measurements 
of  self-  and  latent  hardening  (cf.  [96]),  /?GN  is  the  length  parameter  for  normalization, 

^N  =  |(^SGN-m')|  (86) 


is  the  scalar  magnitude  of  the  projected  lattice  couple  stress  associated  with  the  GND  tensor  on 
slip  system  j  (with  the  dot  product  here  denoting  a  contraction  of  indices  on  the  cotangent  space 
of  bp),  and  SJss  is  the  lattice  stress  associated  with  the  SSDs  on  slip  system  j.  The  strain  measures  of 
SSDs,  f?ss,  are  defined  by  an  additive  decomposition  of  the  statistically  stored  line  densities  on 
each  slip  system,  p‘ss: 


The  corresponding  thermodynamic  force  is  then,  from  (58)4,  (84),  and  (87), 


01 jj 


iSoe  = 

88  de, 


ss 


diA  dess 
dess  dess 


d<A  dess 
dess  dess 


Cifi  yfj\ e'ss  =  c^pJ'  *8sS, 


(88) 


C3t‘hs 
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providing  a  hardening  contribution  in  (85)  proportional  to  the  square  root  of  the  dislocation  line 
density  per  unit  volume,  as  originally  suggested  by  Taylor  [97]  and  since  verified  experimentally 
numerous  times  (cf.  [56,98]),  and  used  in  both  local  crystal  plasticity  theories  (cf.  [99,100])  and 
gradient-based  theories  [34,39,37],  The  GND  (scalar)  couple  stresses  are  found  explicitly  in 
terms  of  the  tensorial  density  of  GNDs  by  inserting  (58)3  into  (86): 


Sgn  = 


rj'x 


Fl 


dAeoiP 


nr 


mx 


—  PgnI1^' 


i-l  ,ps 


F^Al.F^mU  =  ft 


1lta1  .p 


(89) 


mA’ 


withTGN  a  scalar  measure  of  GNDs  on  slip  system  i.  Substituting  Eqs.  (88)  and  (89)  into  Eq.  (85) 
then  gives 


k  —  k0  +  fiJ  (y/fahWJA>w  3“  C^SSj£SS^  •  (90) 

The  contribution  of  GNDs  to  individual  glide  system  resistances  was  previously  suggested  by 
Acharya  and  Bassani  [33,73],  Acharya  and  Beaudoin  [101],  and  Bassani  [74],  for  rate-dependent 
and/or  rate-independent  single  crystal  plasticity,  similarly  to  the  second  term  in  our  Eq.  (90). 

Many  have  used  the  conjugate  stress  to  a  tensorial  GND  measure  as  a  contribution  to  a 
backstress  on  each  slip  system,  rather  than  a  contribution  to  slip  system  resistance  k  as  in  our 
Eq.  (90),  with  a  variety  of  different  ways  proposed  for  projecting  the  tensorial  backstress  onto 
individual  glide  systems  [9,28,32,34,36,37,101].  If  both  a  backstress  and  a  glide  stress  are  desirable 
in  order  to  completely  characterize  the  hardening  behavior  (cf.  [102]),  then  we  can  generalize  the 
flow  rule  in  Eq.  (80)  to 


f  =  isgn  (p‘  -  /) 


(91) 


where  increases  in  the  friction  stress  k‘  >  0  now  manifest  only  from  SSDs: 

k1  =  k0  +  (c3fJ'  ')hlss/ss,  (92) 

and  where  the  backstress  not  necessarily  positive  in  sign,  depends  solely  upon  the  density  of 
GNDs  (cf.  [34]): 

(93) 

In  Eqs.  (92)  and  (93),  h'SSJ  and  h‘GNj  are  slip  system  interaction  matrices,  different  in  value,  but  not 
in  function,  from  those  introduced  already  in  (85).  Different,  and  perhaps  more  inclusive,  ways  of 
expressing  the  contributions  of  SGN  to  hardening  on  each  slip  system  than  the  direct  projection 
method  of  Eq.  (86)  also  merit  further  exploration  (see  e.g.  [28,35,36]). 

While  the  density  of  geometrically  necessary  dislocations  can  be  calculated  directly  from  spa¬ 
tial  gradients  of  the  elastic  lattice  deformation  (Eq.  (49)),  the  densities  of  statistically  stored 
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dislocations  on  each  slip  system  are  modeled  here  as  internal  state  variables  (ISYs)  equipped  with 
separate  equations  dictating  their  evolution.  By  definition,  the  time  rate  of  change  in  dislocation 
density  for  a  given  slip  system,  p'ss,  and  the  rate  of  the  corresponding  dimensionless  lattice  strain 
measure  e'ss  for  that  slip  system  are  related  by 


2£ss  — 


(94) 


A  general  evolution  equation  for  each  dislocation  density  p‘ss  is  now  suggested: 

Pss  =  PkdSish  (4)-  M,  S1)  (J  =  1, . . . ,  n).  (95) 

As  is  clear  from  (95),  the  dislocation  density  rate  on  a  particular  system  i  generally  includes 
contributions  from  the  slip  rate  on  the  y-th  system,  via  the  p'=l  dependence,  as  well  as  the  con¬ 
tributions  of  the  slip  rates  of  other  systems,  through  the  p1^’  dependence.  Likewise,  the  rate  of  p'ss 
may  be  influenced  by  the  current  statistically  stored  densities  on  all  slip  systems  through  the  {SJss} 
dependency  [80,91].  Since  interactions  between  SSDs  and  GNDs  are  not  ruled  out  (e.g.,  trapping 
or  annihilation),  we  also  allow  the  densities  of  GNDs  on  all  slip  systems  to  enter  the  evolution  law 
for  tp  given  system,  through  the  {SJGN}  dependency.  Finally,  the  dependency  of  p'ss  on 
S1  =  F1_1S‘F1_t,  the  pull-back  to  bp  of  the  force  conjugate  to  the  stretch  associated  with  micro¬ 
scopic  heterogeneity  within  the  volume  element,  is  included  since  the  formation  of  intragranular 
cellular  structures  will  likely  spur  dislocation  generation  and  annihilation  at  misoriented  subgrain 
boundaries  [60],  Notice  also  that  F1  will  implicitly  affect  the  density  of  GNDs  through  its  presence 
in  the  multiplicative  decomposition  of  the  total  deformation  gradient,  Eq.  (48).  More  specific  rate 
equations  for  the  SSDs,  with  dislocation  populations  further  divided  into  mobile  and  immobile 
parts  [4,99,100]  and  accounting  for  hardening  and  recovery  on  each  system  at  different  stages  of 
plastic  deformation  (cf.  [34,39])  are  also  envisioned  as  possible  implementations  of  (95). 

With  the  evolution  of  FP  and  the  strain  hardening  prescribed  by  Eqs.  (74)  (95),  we  now  turn 
our  attention  to  the  evolution  of  the  F1  tensor.  The  following  general  evolution  equation,  referred 
to  configuration  bm[,  is  now  proposed  for  the  rate  of  stretching  attributed  to  F1: 

V  =  v'(P,  S',  SGN,  §ss),  (96) 

and  a  similar  dependency  is  proposed  for  the  rate  of  rotation  arising  from  F1,  i.e., 

W1ERRiT  =  Wi(P,S1,SGN,Sss).  (97) 

More  specific  forms  of  (96)  and  (97)  have  not  yet  been  developed.  Elowever,  previous  computa¬ 
tional  micromechanical  solutions  indicate  that  the  magnitude  of  V1  should  remain  small  in 
comparison  to  the  magnitudes  of  both  the  total  applied  strain  and  the  volume-averaged  plastic 
strain  in  deforming  ductile  FCC  single  crystals  [88],  The  normalized  stretch  associated  with  V1  was 
non-negligible,  however,  and  did  attain  a  magnitude  comparable  to  that  of  the  recoverable  elastic 
strain  Cc  in  these  calculations, jis  shown  in  Fig.  6(b).  Guidance  in  formulating  Eq.  (96)  can  also 
stem  from  the  contribution  of  V1  to  the  stored  energy  of  cold  work,  as  demonstrated  in  the  data  of 
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Fig.  6(a)  and  suggested  in  the  Helmholtz  free  energy  function  (68).  We  also  must  note  that  the 
calculations  used  to  generate  data  for  Fig.  6  employed  local  crystal  plasticity  theory  within 
subgrains,  so  the  influence  of  GNDs  as  defined  in  the  present  work  was  neglected  in  that  study 
[88]. 

While  we  currently  lack  any  specific  data  for  the  evolution  of  the  rotation  tensor  R1,  additional 
motivation  for  its  inclusion  in  our  crystal  plasticity  model  is  acquired  by  considering  the  de¬ 
composition  of  the  total  vorticity  in  the  intermediate  configuration  bmt: 


(Dskew  =  W  =  (RcTRe)skew  +  (RcT Ve~ 1 VCRC) skew  +  (V  V 


+  (V1RiRlTV1~1' 


/skew 


/skew 


+wp, 


eWc 


qv'W'v-V, 


(98) 


where  all  skew-symmetrization  is  conducted  with  respect  to  fire  metric  g.  Assume  for  the  moment 
that  plastic  spin  Wp,  spin  due  to  residual  stretching  rate  V  ,  and  total  spin  W  are  prescribed, 
respectively,  via  the  evolution  equations  for  crystalline  slip  (Eqs.  (74)  and  (80)  or  (91)),  the 
evolution  equation  for  strain  due  to  subdivision/heterogeneity  (Eq.  (96)),  and  the  displacement 
boundary  conditions.  Furthermore,  assume  that  second  term  on  the  right-hand  side  of  Eq.  (98), 
the  elastic  spin  due  to  Ve,  which  is  known  through  an  appropriate  rate  form  of  the  constitutive 
relationship  (63)  between  the  applied  stress  and  elastic  strain,  is  negligible  in  comparison  to  the 
other  terms  in  (98)  by  the  typical  assumption  of  small  elastic  strains  in  engineering  metals.  As 
mentioned  previously,  since  V1  is  expected  to  not  exceed  the  same  order  of  magnitude  reached  by 
Ve,  neglecting  the  third  term  in  (98)  may  also  be  a  valid  assumption.  Under  these  assumptions,  the 
skew  elastic  spin  Wc  is  then  found  as 

We  ~  W  —  Wp  -  (V'W'U1)^.  (99) 

Recall  that  Fe  provides  the  current  configuration  orientations  of  the  “average”  slip  plane  normals 
and  slip  directions  in  our  crystal  plasticity  framework  by  virtue  of  Eq.  (75);  i.e.,  texture  evolution 
is  essentially  specified  by  the  time  rate  of  Fe  (cf.  [103]).  With  the  assumption  of  small  recoverable 
elastic  strains,  we  have  Fe  «  Re,  and  the  time  rate  of  elastic  rotation  dictates  texture  evolution. 
Essentially,  we  are  at  liberty  to  completely  control  texture  evolution  through  prescription  of  W1, 
which  in  turn  can  be  thought  of  as  the  contribution  to  the  total  spin  from  intragranular  dislo¬ 
cation  substructure  formation  (e.g.,  grain  subdivision  processes).  For  example,  if  we  specify 
W1  =  (V1_1(W  -  Wp)V1)skew,  then  texture  evolution  (i.e.,  elastic  spin)  will  be  largely  precluded 
d ue  to  substructure  development  within  the  crystalline  volume  element.  On  the  other  hand  if  we 
set  W1  =  0  for  our  evolution  law  (97),  then  crystal  lattice  orientations  will  evolve  as  in  the  classical 
crystal  plasticity  theory.  Numerical  simulations  [50,104]  have  shown  that  the  classical  theory,  with 
W1  =  0,  gives  texture  predictions  that  are  consistently  too-sharp  for  certain  FCC  polycrystals,  at 
least  within  the  framework  of  Taylor’s  [105]  constraints.  So  by  correlating  potential  evolution 
Eqs.  (97)_with  experimental  texture  measurements,  one  could  conceivably  formulate  evolution 
laws  for  W1,  as  speculated  by  Butler  and  McDowell  [50].  Since  grain  subdivision  and  substructure 
development  are  energetically  favorable,  and  hence  expected,  for  some  materials  even  under 
application  of  macroscopically  uniform  boundary  conditions  [7,43,98]  V  and/or  W1  should  attain 
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non-zero  values  during  some  point  in  the  history  of  uniform  deformation  applied  to  such  ma¬ 
terials. 

In  closing  the  discussion  of  evolution  equations,  we  emphasize  that  when  completely  devel¬ 
oping  Eqs.  (95)  (97),  which  we  have  left  quite  general  in  this  work,  one  should  also  consider  the 
restrictions  on  thermodynamic  admissibility  stemming  from  the  complete  dissipation  inequality, 
Eq.  (67).  While  we  have  included  terms  associated  with  dislocation  defect  densities  and  residual 
elasticity  in  the  free  energy  function  (68),  in  most  metals  their  total  cumulative  contribution  to  the 
first  and  second  laws  of  thermodynamics,  (54)  and  (55),  will  remain  relatively  small  in  magnitude 
compared  to  that  of  the  plastic  dissipation  due  to  dislocation  motion  (77),  in  agreement  with 
experimental  measurements  of  the  stored  energy  of  cold  working  [106,107],  This  observation  is 
fully  consistent  with  the  framework  presented  here  (cf.  Eq.  (68));  the  contribution  of  incompat¬ 
ibility-  related  tensors  (i.e.,  F1  and  defect  densities)  to  the  free  energy  is  small  when  averaged  over 
the  entire  volume  of  the  crystal(s).  The  relatively  small  magnitude  of  this  energetic  contribution  is 
scaled  correctly  by  the  appropriate  choice  of  material  parameters  in  (68).  However,  in  our 
opinion,  the  residual  elastic  energy  associated  with  crystal  defects  and  elastoplastic  incompati¬ 
bility  should  not  be  neglected,  as  its  local  density  may  attain  much  greater  magnitudes  in  the 
vicinity  of  local  heterogeneities  such  as  grain  or  phase  boundaries,  and  its  release  may  facilitate 
void  growth,  fracture,  shear  localization,  recrystallization,  and/or  phase  transition  processes  in 
many  metals  [88,108], 


5.  Conclusions 

We  have  used  explicit  volume  averaging  procedures  (Sections  2  and  3)  to  motivate  a  continuum 
formulation  of  gradient  crystal  plasticity  (Section  4).  Notable  features  include 

•  A  three-term  multiplicative  decomposition  for  the  deformation  gradient,  F  =  FeF1FP,  with  Fe 
associated  with  the  average  applied  stress  and  total  lattice  rotation,  F1  accounting  for  the  pos¬ 
sible  development  of  heterogeneous  dislocation  substructures  and  residual  elasticity,  and  Fp 
accounting  for  the  history  of  plastic  deformation  due  to  dislocation  motion. 

•  Prescription  of  a  resolved  Eshelby-type  stress  measure  as  a  driving  force  for  plastic  shearing 
rates,  in  order  to  ensure  a  positive  dissipation  contribution  from  the  plastic  velocity  gradient. 

•  Slip  system-level  strain  hardening  dependent  upon  conjugate  thermodynamic  forces  to  densities 
of  GNDs  and  SSDs. 

•  Allowance  of  F1,  the  kinematic  measure  of  intragranular  heterogeneity  and  possible  subdivi¬ 
sion,  to  influence  the  evolution  of  dislocation  densities,  and  vice-versa. 

•  Allowance  of  a  measure  of  residual  elasticity — the  stretch  associated  with  F1 — to  account  for  a 
fraction  of  the  stored  energy  of  cold  work,  influence  the  effective  elastic  moduli  in  the  unloaded 
configuration,  and  provide  a  residual  free  energy  term  biased  by  the  applied  stress. 

Our  framework  offers  powerful  tools  for  addressing  multiscale  issues  in  plasticity,  most  notably 
providing  a  very  rigorous  (i.e.,  mathematically  precise)  methodology  for  characterizing  kinematics 
and  thermodynamics  with  changes  in  scale  of  observation.  Future  work  should  focus  on  deter¬ 
mination  of  more  precise  evolution  equations  for  F1  and  the  dislocation  densities,  determination 
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of  the  needed  material  parameters  (e.g.,  hardening  coefficients  and  length  scale  constants),  and 
exploration  of  boundary  conditions.  Micromechanical  simulations  of  representative  defect  pop¬ 
ulations  invoking  discrete  dislocation  dynamics  relations  [109-111]  are  foreseen  as  valuable  tools 
for  deducing  evolution  laws  and  model  parameters.  Extensions  to  other  defects  (e.g.,  disclinations, 
point  defects,  and  damage),  temperature  effects,  and  dynamic  conditions  are  also  envisioned. 
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